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The symmetrical restricted Gibbs ensemble (RGE) is a version of the Gibbs ensemble in which par- 
ticles are exchanged between two boxes of fixed equal volumes. It has recently come to prominence 
because - when combined with specialized algorithms - it provides for the study of near-coexistence 
density fluctuations in highly size-asymmetric binary mixtures. Hitherto, however, a detailed frame- 
work for extracting accurate estimates of critical point and coexistence curve parameters from RGE 
density fluctuations has been lacking. Here we address this problem by exploiting an exact link be- 
tween the RGE density fluctuations and those of the grand canonical ensemble. In the sub-critical 
region we propose and test a simple method for obtaining accurate estimates of coexistence densi- 
ties. In the critical region we identify an observable that serves as a finite system size estimator for 
the critical point parameters, and present a finite-size scaling theory that allows extrapolation to 
the thermodynamic limit. 



I. INTRODUCTION 

In a recent paper Liu et al [1J have presented a 
method for studying density fluctuations in highly size- 
asymmetric binary fluids. The method draws on the ge- 
ometrical cluster algorithm (GCA) of Liu and Luijten [2 
which facilitates rejection-free Monte Carlo simulations 
of highly asymmetric mixtures via large-scale collective 
updates which move whole groups of particles (both large 
and small) in a single step. In its basic form, the GCA is 
inherently canonical, i.e. it conserves the total number of 
particles in the system. To enable density fluctuations - 
and thus facilitate the study of phase behaviour - it was 
modified in ref. pQ so that cluster updates transfer parti- 
cles between two subsystems of fixed equal volumes, thus 
realizing complementary density fluctuations within each 
subsystem. Such a scenario constitutes a special case of 
the well known Gibbs ensemble (GE) [3], which has been 
dubbed the restricted Gibbs ensemble (RGE) [H[5] due to 
the absence of volume transfers between the boxes. The 
lack of volume fluctuations in the RGE (which are incom- 
patible with cluster updates) means that the two subsys- 
tems do not remain at equal pressure - a fact which com- 
plicates the determination of phase coexistence points. 

In this paper we address the issue of how to determine 
coexistence curve and critical point parameters within 
the RGE. For simplicity we perform the analysis for a 
single component fluid, but we expect the results to be 
directly applicable to the asymmetrical binary case in 
which the GCA+RGE method operates [2T]. Our refer- 
ence point is the well studied grand canonical ensemble 
(GCE) for which powerful techniques exist for determin- 
ing critical points and coexistence curve properties [BJ. 
In order to develop methods of comparable utility, we 
appeal to an exact link between the density fluctuations 
in the RGE and those of the GCE. This link is exploited 
to predict the properties of simulations performed in the 



RGE simulations on the basis of data accumulated from 
GCE simulations. A detailed analysis shows, in partic- 
ular, how one can locate criticality and the coexistence 
binodal within the RGE. 



II. LINKING THE RESTRICTED GIBBS AND 
GRAND CANONICAL ENSEMBLES 

The symmetrical restricted Gibbs ensemble comprises 
two subsystems (or boxes) of fixed equal volumes V± — 
V 2 = V = L d (with d = 3 in the simulations to be con- 
sidered below) containing respective particle numbers iVi 
and N 2 at some imposed temperature T. The boxes can 
exchange particles subject to the constraint that the total 
particle number 



Ni + N 2 = N 



(1) 



is fixed. The probability of finding such a restricted 
Gibbs ensemble system in a given microstate having par- 
ticle numbers iVi^ and coordinates {q}i,2 is given by 

P"iN,,.{,}„\N»V,T) = ^-** + ™ . (2) 

where /3 = 1/ksT is the inverse temperature, while E\ 2 
is the configurational energy of the respective boxes. 
From this it follows that the probability of finding the 
system with given N\ (and hence prescribed N 2 ) is 

P RG (N\N V T)- Z(NuV,T)Z(N 2 ,V,T) 

P (N^No^T)- zRG ( No}VjT) ( 3 ) 

where Z(N, V, T) is the canonical partition function. 

Now to forge the link to the grand canonical ensemble, 
we note that 

*W = »B, (4) 



2 



where P(N\p,V,T) is the grand canonical probability 
distribution at chemical potential p. Inserting into eq. ^ 
yields 

(N 1 \N ,V,T) = wP(N 1 \ t x,V,T)P(N 2 \fx,V,T) (5) 



RG 



P 



where the transformation is defined for all N\ < No and 
the normalization constant, w, is given by 
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[Z{p,V,T)} 2 exy{-{3pN ) 
Z« G (N ,V,T) 



(6) 



Eq. (|5| is formally exact and provides the necessary 
link between the grand canonical and symmetrical re- 
stricted Gibbs ensembles [55]. An interesting property 
of P RG is its independence of the choice of the chemical 
potential appearing on the right hand side of Eq. ([5])- a 
feature that has its origin in the complementarity of the 
fluctuations of P{Ni) and P(N — N\) |2"g] . 

Of course, it is more conventional to work with the 
number density rather than the particle number, in which 
case the transformation takes the form 



Pl°{pW T) = P L (p\T)P L (2 Po - p\T) , 



(7) 



where p — N\/V is the single box density, Pl(p\T) is 
its GCE distribution at temperature T, p — N n /2V is 
the fixed overall density in the RGE, and = means to 
within an unspecified normalization constant. For clarity 
we have suppressed reference to the arbitrary chemical 
potential p and the constant volume V, but have added 
a subscript L to emphasize that we nevertheless remain 
concerned with systems of finite size. 

A central feature of P^ G (p) is its symmetry: the trans- 
formation ([7 1 involves the product of Pl(p) and a re- 
flected and shifted (by 2p ) version of itself. Quite gen- 
erally [and irrespective of the symmetry of Pl(p)], this 
implies that P RG (p) is symmetric (i.e. even) with a mean 
value of po. This symmetry (and the p independence of 
P RG (p) that it implies) precludes deployment of the stan- 
dard tools for locating phase coexistence and criticality 
that have been developed in the context of the GCE [6J . 
Instead bespoke analysis techniques are called for. 

In this paper we develop such techniques. Perhaps sur- 
prisingly, our strategy does not involve performing actual 
RGE simulations. Instead we apply the exact transfor- 
mation Eq. ([7]) to independently obtained GCE density 
distributions for the Lennard- Jones fluid in order to in- 
fer the properties of the resulting RGE density distribu- 
tion. Our rationale for so doing is that it is well under- 
stood how one determines coexistence and critical point 
properties from GCE density distributions. By starting 
from such, we have a known baseline with which to com- 
pare the results of our analysis of the RGE. Accordingly, 
within this approach, GCE distributions serve as input 
to Eq. ([7]) and pa simply becomes a parameter of the 
transformation. As it is varied, for some prescribed T, a 
spectrum of symmetrical distributions P RG (p\po) is gen- 
erated, each centred on the respective value of po- In sec- 



points in the neighbourhood of liquid-vapor coexistence 
and present a method for extracting coexistence densi- 
ties. Thereafter, in section|V] we consider state points in 
the vicinity of the liquid- vapor critical point, and perform 
a finite-size scaling analysis of P^ G (p) which facilitates 
accurate estimates of critical point parameters. To be- 
gin with, however, we shall briefly introduce our model 
system. 



III. FLUID MODEL 

The GCE density distributions for use as input to the 
transformation ^ were obtained via Monte Carlo simu- 
lations [7] of a 3d fluid interacting via a Lennard- Jones 
(LJ) potential: 



cp( r )=Ae[{a/r) 12 -{a/rf 



(8) 



Here r is the particle separation, a sets the length scale 
and e is the well depth. The potential was truncated at 
r c = 2.5a and no correction applied, thus allowing com- 
parison with the results of a previous study performed 
under the same conditions [BJ. As is customary, we shall 
refer to the reduced temperature T = l/(/3e) rather than 
the well depth in the simulation results to be described 
below. Furthermore, we shall normally express tempera- 
ture as a multiple of the critical value T c = 1.1878(2) (as 
determined in Sec. [V}. Further details of the simulation 
methodology can be found in ref. [BJ. 

For the studies of the subcritical region described in 
Sec. |IV| a single periodic system of size L — 10a was em- 
ployed. The finite-size scaling investigation of the critical 
region, described in Sec. |V| used seven periodic systems 
of linear extent L = 10a, 12.5cr, 15c, 17.5a, 20a, 22.5a, 
25a. We set a as the unit of length below. 



IV. SUBCRITICAL COEXISTENCE REGION: 
THE INTERSECTION METHOD 

We commence by addressing the behaviour of the RGE 
density distributions in the subcritical coexistence region 
of our model fluid. As described above, our approach will 
be to glean information about the RGE not via explicit 
RGE simulations, but by applying the exact transfor- 
mation Eq. ([7]) to GCE simulation data. Suitable esti- 
mates of the GCE density distribution P^ (p) have been 
obtained in a previous study of the Lennard- Jones fluid 
by one of the authors [BJ. Examples for temperatures 
T = 0.842T C and T = 0.69T C are shown in Fig. [TJ in 
which the chemical potential has been tuned so that the 
peaks have equal areas - the criterion for coexistence [H] . 

The two peaks in the coexistence form of Pl (p) derive 
from fluctuations around the respective densities p g and 
Pi of the gas and liquid phases. These densities are found 
as an average over the appropriate peak in Pl{p) via 



tion IV we consider the nature of this spectrum for state 



pP L (p\T)dp. 



(9) 
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FIG. 1: The form of the GCE density distribution P L (p) for 
temperatures T — 0.842T C and T = 0.69T C , expressed on a 
logarithmic scale. 



where p m i n and p max are chosen suitably to fully encom- 
pass the appropriate peak of P^ G (p), ie p m i n = p g .i — S, 
Pmax = Pg.i + S. For our purposes, a useful quantity will 
prove to be the average of the coexistence densities at a 
given temperature, known as the coexistence diameter : 



Pd {T) = \p g {T) + Pl (T)]/2. 



(10) 



A pertinent feature of Fig. |T]is the shape of the distri- 
bution in the central region of density, between the peaks. 
Here one sees that the distribution flattens markedly, par- 
ticularly so for the lower temperature. This flattening 
corresponds to the appearance of well defined interfaces 
between the coexisting phases. For sufficiently low T 
or sufficiently large system size L, the central portion 
becomes completely flat because the interfaces decouple 
from one another, there being no free energy penalty for 
one phase to change its fractional volume with respect to 
the other [9]. 

Let us now consider typical forms of P RG (p\po) and 
how they derive from the underlying GCE density distri- 
bution. Examples are shown in Fig. [2] for widely sep- 
arated values of po- Specifically, the figure plots the 
form of Pf G (p) for T = 0.842T C at three values of 
po = 0.07, 0.30, 0.6. These show that for the smallest and 
highest values of po, P RG (p) is single peaked, while for 
the intermediate value of po, P^ G (p) is double peaked, 
though in this latter case neither peak coincides with ei- 
ther the gas or liquid densities. 

A fuller picture of the po dependence of P RG (p) 
emerges by considering the peak densities as a function 
of po- Let us denote the densities of the low and high 
density peaks of P^ G in the double-peak regime as p_ 
and p+ respectively, these being obtained as the average 
RGE peak density in a manner analogous to Eq. A 
plot of p_ and p + versus po is shown in Fig. [3] for the 
temperature T = 0.842T C . From this figure one sees that 
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FIG. 2: The form of Pl° (p\po) (expressed on a log scale) for 
the LJ fluid at T = 0.842T C for (a) p = 0.07; (b) p = 0.30; 
and (c) po = 0.6. In each instance the dashed curve is the 
coexistence form of the GCE distribution, Pl(p); the dotted 
curve is the distribution, P^(2p — p), and the solid curve is 
the resulting RGE distribution, Pl°{p) obtained via Eq. d7|. 
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FIG. 3: The peak densities of Pl G {p) as a function of po 
for the Lennard Jones fluid with T = 0.842T C . The dotted 
lines on the ordinate at low, intermediate and high densi- 
ties correspond to the values of p g , pd and pi respectively, as 
determined independently from the underlying GCE data as 
described in the text. The value of pd is also marked by a 
vertical dot-dashed line on the abscissa. 



as a function of po, the RGE distribution is single peaked 
at low density, splits into two peaks at intermediate po, 
but that these peaks merge back into a single peak at 
high values of pq. Furthermore in the two-peak regime, 
p_ and p + vary considerably, though we note that they 
remain subject to the constraint 



Pa = (p-+p+)/2 



(11) 



which is mandated by the symmetry of P RG (p), together 
with the fact that its average is strictly equal to po. 
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TABLE I: The combinations of p (rows) and 2po — p (columns) 
corresponding to p g , pd and pi . The entries in the table give 
the resulting po. 



Now, it transpires that there are certain special values 
of p for which p_ and p + both coincide with one of p g , 
Pd and pi. To substantiate this, consider the derivative 
of P HG {p). From Eq. (T7|). This derivative certainly van- 
ishes when P' L (p) — P'ij2po — p) = 0. Assuming Pl{p) to 
be symmetric, this happens when both p and 2p — p are 
one of pg, pd, pi- The various possibilities are shown in 
table |IV| where rows are labelled by p and columns are la- 
belled by 2/?o — P, and the entries are the resulting values 
of Pq. Restricting attention to the double peaked regime 
of P RG , one finds from the table the following three re- 



lationships between the peak positions in the RGE and 
GCE ensembles: 



P- = Pg 
P+ = Pi 



when po = p d 



(12) 



and 



P- 
P+ 



P- 
P+ 



Pg 
Pd 



Pd 
Pi 



when pq = (p g + p d )/2 , (13) 



when po = (p t + p d )/2 . (14) 



To demonstrate that these predictions are consistent 
with the data of Fig. [3] we have marked on the ordinate 
the coexistence densities p g , pi and the coexistence diam- 
eter pd, while the abscissa is marked with the densities 
Po = (p g +Pd)/2,po = pd, andpo = (pi+p d )/2. Inspection 
of the figure confirms the above results. Additional in- 
sight and further confirmation is furnished by comparing 
the forms of P RG (p) with the underlying GCE distribu- 
tions for pq = pd and for po = (p g + pd) /2, as shown in 
Fig.0 

Eq. ( fl2| implies that if one can determine the coex- 
istence diameter from a RGE simulation, then one can 
simply read off the coexistence densities from the corre- 
sponding peaks in P RG (p\pd). In practice, one can esti- 
mate pd by requiring consistency between Eq. ( |12| and 
Eq. ( fl3| and/or Eq. ( fl4| . The consistency condition is 
readily implemented graphically as shown in Fig. |5j one 
simply plots the measured values of p± versus po, to- 
gether with the transformed points p± versus 2po — p± . 
The two sets of data intersect at p = Pd and the coex- 
istence densities can be read off from the position of the 
RGE peaks at this value of pq [24]. The implementation 
of this "intersection method" is shown in Fig. [5| 

In order to gauge the accuracy of our intersection 



method, we tabulate in Tab. IV the results of applying 
it at a number of subcritical temperatures spanning the 
range 0.7 < T/T c < 0.94. Clearly there is excellent agree- 
ment between the coexistence densities obtained from the 
intersection method and those obtained directly from the 
underlying GCE coexistence density distribution. Such 
discrepancies as are evident are largest at temperatures 
close to criticality where the positions of peaks in a den- 
sity distribution anyway serve as a poor measure of co- 
existence densities due to finite-size effects [6]. 

Finally in this section, we note that in the regime of 
system size and temperature for which a minimum occurs 
between the peaks of Pl (p) (rather than a flat portion) , 
the rigorous validity of Eqs. ( fl3) and ( J14] ) rests upon the 
coincidence of the coexistence diameter with the density 
of this minimum, i.e. on the symmetry of Pl{p)- In gen- 
eral we observe very close, but not perfect, agreement be- 
tween these two densities. However, in practice it seems 
that any deviation of the two quantities will impact lit- 
tle on the effectiveness of the intersection method. This 
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RGE 




T/T c 


Pd 


Pa 


Pi 


Pd 


Pa 


Pi 


0.937 


0.3320(16) 


0.1095(12) 


0.554(4) 


0.3352(14) 


0.1096(8) 


0.5608(22) 


0.865 


0.3473(24) 


0.0607(8) 


0.6342(44) 


0.3482(20) 


0.0600(10) 


0.6363(44) 


0.844 


0.3517(22) 


0.05079(74) 


0.6526(44) 


0.3513(24) 


0.04998(66) 


0.6527(42) 


0.833 


0.3538(26) 


0.0467(8) 


0.6611(50) 


0.3537(22) 


0.04598(50) 


0.6615(42) 


0.813 


0.3582(12) 


0.0396(12) 


0.6769(26) 


0.3583(16) 


0.0390(4) 


0.6777(32) 


0.796 


0.3619(30) 


0.0342(2) 


0.6896(60) 


0.3616(26) 


0.03378(30) 


0.6895(54) 


0.695 


0.3847(12) 


0.01360(8) 


0.756(2) 


0.3854(14) 0.013698(90) 0.7572(26) 



TABLE II: Estimates of the coexistence densities p g and pi, together with the coexistence diameter pd for the 3D LJ fluid 
described in Sec. |III| at a selection of sub-critical temperatures. The GCE estimates of p 9 ,; derive from the densities of the 
peaks in the coexistence form of Pl(p); Pd follows as their average. The RGE estimates for pd derive from Pl°(p) vm the 
intersection method described in the text, while p g> i follows from the corresponding peak densities in Pl G (p\po = Pd)- Error 
bars are conservative and are calculated via a blocking analysis of the underlying GCE data. 
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FIG. 5: Illustration of the intersection method. Circles show 
the RGE peak densities p_ and p+ of Fig. [3] as a function 
of po- Also shown are the transformed points, 2po — p±, 




(squares). The two data sets intersect at the diameter (verti- 
cal dot-dashed line). 



dependent of po for a range of values of po w p d . Under 
these circumstances the coexistence densities can simply 
be read off as the values of these constant densities. 



V. NEAR-CRITICAL FINITE-SIZE SCALING 



FIG. 4: (a). The solid curve shows the form of Ff a (p) for 
the case po ~ Pd, for which p_ = p g and p + = pi. Also 
shown (dashed curve) is the underlying GCE coexistence in- 
put distribution Pl(p) and the reflected and shifted distribu- 
tion Pt(2po — p) (dotted line) that feature in the transfor- 
mation Eq. 171. (b). The corresponding data for the case 
Po = {pg + Pd)/2, for which p_ = p g and p + = p d . In both 
parts the coexistence diameter is represented by a vertical 
dot-dashed line. 



is because near its minimum, Pl{p) is generally much 
more slowly varying than near its peaks. To elaborate, 
consider the case po = (p g + Pd)/% shown in Fig. Qb). 
From the figure one sees that the lower density RGE 
peak is formed by multiplying the peak in Pl(p) by the 
minimum in P^(2po — p). As the system size increases, 
the peak becomes sharper, the minimum becomes flat- 
ter, and hence any discrepancy between the density of 
the minimum and that of the coexistence diameter will 
have increasingly less influence on the location of the low 
density peak of Pf; (p). The same argument holds for 
the high density peak of Pl G {p) at Po = (Pg + Pd)/^- Of 
course, in the limit in which the central range of Pl(p) 
becomes flat, the values of p_ (and p + ) will become in- 



A. Symmetric fluids 

The near-critical finite-size scaling properties of the 
RGE have been considered in two previous studies. The 
first by Mon and Binder [3] focused attention on deter- 
mining the critical temperature for a lattice gas model, 
the critical density being known a-priori by virtue of 
particle-hole symmetry. Whilst successful with respect 
to its aims, this study provided no insight into how to de- 
termine the critical parameters for off-lattice fluids. The 
second study, by Bruce [5], provided elucidation of the 
scaling form of the critical point density distribution (and 
its relationship to the critical Ising magnetisation distri- 
bution), but here too no method was given for locating 
criticality in realistic fluids within the RGE. 

Here we build on these studies by providing a sim- 
ple prescription for locating critical points for off-lattice 
fluids. To begin with, however, we shall follow Mon and 
Binder's lead by seeking inspiration from a model exhibit- 
ing particle-hole symmetry -namely the Ising lattice-gas 
model - and consider the fluctuations in the density p. 
We shall first discuss the critical limit before moving on 
to the near-critical region. 
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1. The critical limit 

In the critical limit t = (T - T c ) /T c 0, h = (/x - 

Mc)/^c — > 0, and for sufficiently large L, Pl{p) adopts 
the well known universal scaling form [TO] 



P L (p\p, c ,T c ) = a L p*(a L (p - p c )) 



(15) 



Here ai, = a\L^^ v , with a\ a non-universal scale factor, 
while (3 ~ 0.326 and v ss 0.630 are the usual critical expo- 
nents for the order parameter and the correlation length 
respectively [IT], p* is a universal scaling function, which 
is symmetric (even) with respect to the mean density p c . 
Its form is characteristic of the Ising universality class, 
and is well established from previous simulations studies 
of the Ising model in both d — 2 and d = 3 [TU H3] • 

Let us now consider the RGE density distribution 
P RG (p\po,T c ), and make the choice po — p c . Then, from 
the symmetry of the critical grand canonical density dis- 
tribution, we have that Pr,(p\p c , T c ) = Pj,(2pQ — p\p c , T c ). 
Utilizing Eqs. and (15 1 now shows that 



P* G {p\p Cl T c ) = [p*{a L {p- Pc ))f 



(16) 



up to an unspecified normalization factor. Thus the crit- 
ical point RGE distribution matches the square of the 
known universal fixed point form of the Ising magnetisa- 
tion distribution, as also previously pointed out by Bruce 
0. 



2. Deviations from criticality 

At the critical temperature T = T c , the form of P RG 
for values of po m t ne vicinity of p c can be related to 
that for po = Pc by expanding in powers of a reduced 
density go = po — p c - To facilitate ready comparison of 
the various members of the resulting spectrum, it is fur- 
ther convenient to work with a version of P RG (p) having 
zero mean (recall that po not only controls the form of 
this distribution, but also sets its mean). To this end we 
set x = p — p c and define a shifted GCE density distri- 
bution 



V L (x\n c ,T c ) = P L (x + p c \fi c ,T c ) , (17) 



so that 



P HG (x\p ,T c ) = V L (x + Qo\pl c ,T c )V l {-x + Q Q \y. c ,T c ) 



Expanding in the small parameter one then finds 
quite generally that 



P, 



* G (x\p ,T c ) = V L (x)V L (-x) 



+ Qo [P L (x)V' L (-x) + V L {-x)V' L {x)} 
+ \qI[Pl{x)V'1{-x) + 2V' l {x)V' l {-x) 

+ V L {-x)V'[{x)] 
+ \qI[ZV' l {x)VI{-x) + W' l {-x)V'1{x) 

+ V L {x)V'H{-x)+V l l!(x)V L (-x)] 



24 



{8VUx)V'l(-x) 



(19) 



where we have abbreviated Vl(x) = Vl(x\p c ,T c ). 

For models exhibiting particle-hole symmetry, the 
terms in the expansion ( 19 ) that involve odd powers of 
go drop out |25j and one has 



Pi 



fJ G (x\p ,T c ) 



V l {x)Vl{-x) 

qI [v L {x)v'l{x)~{V L {x)f}+o[ e t\. 

(20) 



Thus in the symmetrical case, the leading po dependence 
of P RG (x) comes from the p§ term. An analogous ex- 
pansion of P RG (x) in terms of the reduced temperature 
t shows that, by contrast, the leading temperature de- 
pendence arises from terms linear in t. 

Together these findings lead us to propose the following 
finite-size scaling ansatz for P RG (x\po,T), expressed in 
terms of the reduced variables go and t: 



P, 



L°{x\go,t) 



b L p(b L x, b 2 tl>l\ b 3 g 2 L 2f) /», b,L- 6 ^) 



(21) 

Here b L = b l L fi / u , while b u b 2 MM 

are non-universal 

scale factors, p is a universal scaling function which is 
symmetric in x for all values of its other arguments. The 
final term, with associated exponent 8 « 0.52 [T3] charac- 
terizes the leading symmetrical corrections to scaling and 
accounts for deviations from the large L scaling limit. 

Now, a useful dimcnsionless quantity that character- 
izes the shape of a distribution is the fourth order cumu- 
lant ratio, a common variant of which is defined as 



Q = 



(22) 



The scaling of the cumulant ratio for P RG (x\g , t) follows 
from Eq. pll) as 



Q L (p , t) = q(c 2 tL^ v , c 3 g 2 L 2 P/», c^- e l v ) , (23) 

with q a universal function. Accordingly in the vicinity 
of the critical point the cumulant behaves as 



Ql(Po^) 



1 



(18) 



+ cl i L- 6 '» + 0{glt 2 ) 



(24) 
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where q* = g(0,0,0) and c' 2 ,c' 3 and c 4 are system-specific 
constants. 

For a given L and a fixed choice of T, this implies that 
the variation of Ql with po is parabolic: 



Xp ,t = 0)=q* 



const. + c' 3 L 20/ "qI 



(25) 



This latter result can be checked directly for the case T = 
T c because the universal critical point form p*{ixl{p~ Pc)) 
has been parameterized on the basis of accurate measure- 
ments of the magnetisation distribution of the critical 3d 
Ising model [13]. Fig. |6ja) shows a selection from the 
spectrum of transformed distributions P^ G (x\g ,t = 0) 
for a wide range of values of go- For small values of go, 
the behaviour of Ql(po) is well described by a parabolic 
form with q* = 0.711901, as confirmed by the fit shown 
inFig.gb). 

Finally in this section, it will prove useful to define 
a special locus in (go,t) space, namely that for which 
Ql matches q* , which we dub the "iso-cr*" line. Let us 
initially disregard corrections to scaling (i.e. set C4 = 0). 
Then for points on this line sufficiently close to criticality, 
it follows immediately from Eq. ( 24 1 that 



2 



-(4/4)^-^, 



(26) 



which is a parabola in the space of po — T whose maxi- 
mum coincides with the critical point. Now, the presence 
of corrections to scaling means that at the critical point 
(t = go = 0), Ql will differ from its limiting value q* 
by an amount c' 4 L~ e / v . Effectively, therefore, in follow- 
ing the iso-q* parabola to its maximum for a finite-sized 
system, we miss the critical point by a small deviation. 
Empirically, one finds that c' 2 and c' 3 are negative [cf. 
the effect of finite go in Fig. [(^a)]. Since for Ising-like 
systems one finds that c 4 is negative, i.e. corrections to 
scaling decrease Q l below q* at criticality [BJ [JJ] , it falls 
to negative values of t to raise Ql back up to q*. This 
implies that the maximum of the iso-q* line occurs at 
a temperature T C (L) which differs from the true critical 
temperature according to: 



T c - T C {L) cc L 



-(e+i)/u 



(27) 



Values of T C (L) obtained from the maximum of the iso- 
q* parabola for a number of system sizes can thus be 
extrapolated to yield an estimate of T c . 



B. Asymmetric fluids 




0.715 




FIG. 6: (a) Selections from the universal spectrum of trans- 
formed distributions generated atT = T c via Eqs. H15|18| | us- 
ing the estimate for p* (x) parameterized in ref . [13] . All distri- 
butions are normalised. The scaling variable is x = 61 x, 
with the non-universal constant bi chosen so that for distri- 
bution for go — has unit variance, (b) The corresponding 
behaviour of the cumulant Q(qq) together with a parabolic 
fit. The scaling variable here and in the legend of (a) is 
£0 = ^/csL 0/v g o . 



variables. The magnitude of these corrections (relative to 
that of the limiting form p*) decays with L [T5J [TH1 ITT] . 

Explicitly, one can confirm by a detailed finite-size 
analysis [TJ] the form conjectured by [TJJ, as an expan- 
sion in inverse powers of L and up to quadratic order in 



In realistic fluids, the critical point density distribu- 
tion in the GCE is genuinely symmetric only in the limit 
L — > 00 . For finite-sized systems at criticality, corrections 
are expected from the fact that the Ising-like scaling fields 
are not identical to t and h, but are instead given by lin- 
ear combinations of these two variables ( "field mixing" ) 
and also the pressure ("pressure mixing"), with in addi- 
tion quadratic and higher-order corrections in the same 



a^PLiplPcT) = p* +jL-^» Spm + a 4 L- e ^r cs (28) 
+ lL^- a -^ v s iTa + a 2 L 1 / v tr t 



Here p* 



r t, ^t 2 ,fm are even functions, while s pm , Sf m 



and fin are odd; all are evaluated at the scaled density 
cll(p— Pc)- All terms on the r.h.s. except for p* have non- 
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universal prefactors which we have written as j, k, I, ai 
and (14. The subscripts indicate where the various contri- 
butions come from, according to pm: pressure mixing, cs: 
corrections to scaling, fm: field mixing, t: temperature 
shift, t 2 : second order in temperature shift. 

From Pl(p\[i c) T) as given above one can find 
P^ G (x\po, T), and then take moments of x to obtain the 
cumulant ratio Q. This basically inherits the structure 
of the terms in Pl(p\p c ,T), except for the fact that the 
asymmetry of s pm and Sf m cancels linear contributions in 
the mixing coefficients: 

Q L (po,t) = q^l + d.tL^ + cylL^ + c^L- 6 /" (29) 
+ ...fL- 2 l 3 /» + ...jkL (2 - a - 2 V/»t 

+ ...jg L- 2 ^" + ...kg L( 2 - a - 2 ^} 

The dots (...) indicate factors of order unity that are 
determined by the various universal functions r and s. 
It is then easy to see, by setting dQhjdgo = 0, that the 
extremum of any line of constant Ql in the (go,t) plane 
occurs at 



QO 



.jL- 2 ^" + ...kL^- a - 2 ^^t 



(30) 



rather than at go = as in our previous analysis with- 
out field and pressure mixing. Inserting back into (29) 



then shows that the value of Ql on any iso-Q^ line has, 
as a function of the reduced temperature t at the maxi- 



mum, the same form as the r.h.s. of (29) but with the go 



dependent terms omitted and the dots now representing 
different prefactors of order unity. By setting Ql = q* 
one concludes that the maximum of the iso-g* curve, in 
particular, lies at the reduced temperature obeying 

= c! 2 tL 1 ' 1 ' + J i L- e ' v + ...fL-W" (31) 
+ . . . 3 kL (2 ~ a ~ 2 ^l v t +... Jfe 2 £ 2 ( 2 -«-/»)/"i a 

The corresponding value of go at the maximum follows 
by inserting the solution for t into ( 30 1 . 



We have kept track of pressure mixing effects via the 
coefficient j above, but in non-ionic fluids these are gen- 
erally weak, i.e. j is very small [17 . To see how mix- 
ing effects enter, let us then only keep k for now. If k 
is small, the condition (31) will give the l ead ing scaling 
t ~ L- (e+1)/v as before. The last term in ^ is then of 
order ] t 2 i, 2 ( 1 - a -0- e )/ L ' anc ] becomes comparable to the 
first two terms ~ L~ 6 I V when 



L 



-u/{l-a-P-B/2) 



(32) 



The exponent on the r.h.s. is around -2.07, so if k is 
small enough then this crossover lengthscale above which 
field mixing becomes relevant may be too large for us to 
access. Our numerical results below suggest that this is 
the case for the specific system studied here. 

In Eq. ( 30 1 , the pressure mixing term is generally 



kL^ 2 ~ a ~ 2 ^/ u t. Even if mixing is too weak to be ob- 
served in the L-dependence of t, it will therefore be evi- 
dent in an L-dependence of go- In the specific case where 
t ~ [ J -( e + 1 )/ 1 ' t this dependence would be 



go ~ fc^i—V-')/* (33) 
with exponent « —0.45. 

C. Computational method 

Following the iso-g* parabola to its maximum provides 
a convenient route to the effective critical parameters. In 
practice, this is achieved via the following computational 
prescription. 

1. For a given system size L, perform a RGE sim- 
ulation for some value of po and T to obtain 
P?°(P\P0,T). 

2. Deploy histogram reweighting [18 to locate the 
temperature for which Ql(T\po) = q* = 0.711901. 

3. Repeat for a range of po values to yield the iso-g* 
curve in p a — T space. 

4. A parabolic fit to the estimates for the \so-q* curve 
will pinpoint its maximum or - should it prove nec- 
essary - guide the choice of parameters for further 
simulations nearer to criticality. 

5. The coordinates of the maximum of the iso-g* 
parabola serve as estimates for p c (L) and T C (L). 
Results from a range of system sizes can be ex- 
trapolated to the thermodynamic limit according 
to Eqns. p3l and p7l. 



With regard to this procedure, we make the following 
observations. Firstly, our approach explicitly assumes 
that the universality class and the associated value of 
q* are known a-priori. In most cases of current inter- 
est, fluid criticality appears to be Ising-like, and hence 
the above procedure should provide a precise and effi- 
cient procedure for estimating T c and p c . However, in 
cases where the universality class is in doubt, or one sim- 
ply wishes to perform a consistency check on one's mea- 
surements, a procedure can be implemented analogous 
to that proposed by Kim and Fisher [2U] for the GCE 
density distribution. By scanning po at some prescribed 
T, one can locate the maximum of the parabola Q(po\T) 
(cf. Eq. (25)). Repeating for a range of T (via histogram 



subleading because of its L-dependence, so that go 



reweighting) allows the determination of the line of Q 
maxima in po — T space. Estimates of Q along this line, 
and for a range of system sizes, should exhibit "cumu- 
lant crossing" [TU] at the estimated critical point, and 
for a value of Q = q* characteristic of the appropriate 
universality class. Correspondingly the finite size forms 
of P^ G (x) should all collapse onto a universal scaling 
function. 
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-0.0025 



-0.01 

P -P C (L) 



0.0 1 



FIG. 7: Measured iso-g* curves for the LJ fluid for system 
sizes L = 10, 15, 20. Uncertainties are considerably smaller 
than the symbol sizes. Also shown (curves) are parabolic fits; 
maxima provide estimates of the critical temperature T C (L) 
and density p c (L). 



D. Results 

In order to investigate the finite-size scaling properties 
of the RGE, we have generated near-critical forms of P^ G 
from GCE density distributions via the transformation 
Eq. 0. 

Measurements of the \so-q* curve for the system sizes 
L = 10, 15, 20 are shown in Fig. [7ja) together with 
a parabolic fit. Clearly the data are well described by 
a parabolic form, but the positions of the maxima vary 
considerably. 

In order to expose the scaling behaviour of the iso- 
q* curves we have shifted the data so that the maxima 
coincide with the origin, as shown in Fig. [8ja) . We then 
scale the temperature axis by L^~ 2l3 ^ v t, as predicted by 
Eq. (26 1 to obtain the results shown in Fig. [8Kb) . This 



o.ooo 



indeed shows an excellent data collapse. 

The coordinates of the iso-q* maxima serve as esti- 
mates for the effective critical parameters T c (L),p c (L). 
The finite-size dependence of these estimates are plotted 
in Figs. |9ja) and (b) respectively in terms of the antici- 
pated scaling variables (cf. Eqs. (27 1 and (33)). As antic- 
ipated above, data for T C (L) are most consistent with the 
scaling t ~ L~( e+1 '' 1/ , suggesting that field mixing effects 
on this quantity are weak, at least for the system sizes 
we have studied. An extrapolation of the data for T C (L) 
(Fig.ga)) yields T c = 1.1878(2). As regards p c (L), the 
range of accessible system sizes and the rather large error 
bars unfortunately prevent us from establishing convinc- 
ingly the predicted L-dependence (33). The error bars 



were determined from a blocking analysis and one notes 
that those for T C (L) are relatively much smaller than 
those on p c (L). The reason for this are twofold. Firstly 
the iso-g* curves have quite small curvature (cf. the 




P -Pc( L > 



FIG. 8: (a) The iso-g* curves of fig. [7] shifted so that their 
maxima lie at a common origin. Statistical errors are smaller 
than the symbol sizes; lines are guides to the eye. (b) The 
same data plotted in terms of the scaling variable L^~ 2 ^' u t, 
with (3 = 0.326, v = 0.63. 



scales of Fig. [7]) so that any uncertainty in T translates 
to a larger uncertainty in p. A further problem is that for 
the smaller system sizes the resolution of the density scale 
Ap = l/V limits the precision with which the maximum 
can be determined. For instance for L = 10, the den- 
sity resolution is Ap — 0.001 which is large on the scale 
of the finite size shifts we observe overall. If we restrict 
our extrapolation to a fit based only on the five largest 
system sizes shown in Fig. (9jb), we find p c = 0.3210(3). 
The present estimates for T c and p c are to be compared 
with previous older estimates (based on a more limited 
range of system sizes) of T c = 1.1876(3), p c = 0.3197(4). 
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(b) 



P C (D 
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FIG. 9: (a) Measured estimates of T C (L) from the maxima 
of the iso-q* curve. The data is plotted against the scaling 
variable L ( - e+1 ^ v and a straight line fit shows the quality of 
the scaling. Statistical errors do not exceed the symbol sizes, 
(b) Measured estimates of p c {L) from the maxima of the iso- 
q* curves, plotted against the anticipated scaling variable, 
Eq. pi 



VI. SUMMARY AND CONCLUSIONS 

The aim of the work described in this paper was to 
lay firm foundations for future investigations of phase 
behaviour within the Restricted Gibbs ensemble. To this 
end we have developed new approaches for determining 
both the coexistence and critical point parameters within 
this ensemble. 

In the subcritical regime, we have proposed and tested 
an "intersection method" for estimating the coexistence 
densities. This involves measurements of the RGE peak 
densities as a function of the overall system density po. 
Comparison of the resulting data with a transformed ver- 
sion of itself yields an intersection at a density that serves 
as an estimate of the coexistence diameter pd- The coex- 



istence densities are then simply read off from the RGE 
peak densities for po = pd- In tests, the results agreed to 
high precision with independently determined estimates 
of the coexistence properties. 

In the near-critical region, we have described (and ex- 
tended) a finite-size scaling method (elements of which 
were originally outlined in ref. [I) for obtaining accu- 
rate estimates of fluid critical point parameters within 
the RGE. The strength of this method is that it provides 
estimates of effective (finite-size) critical parameters via 
a simple and convenient parabolic fit to measurements 
of a cumulant ratio. These estimates can be extrapo- 
lated to the thermodynamic limit using appropriate sear- 
ings forms which we have presented both for the case of 
particle-hole asymmetry and for asymmetric fluids. We 
expect the method to be useful when one has prior knowl- 
edge of the universality class and the associated value of 
q*. 

In addition to its utility for extracting critical point 
parameters from RGE simulations of fluids, we remark 
that our analysis provides an alternative to traditional 
routes for obtaining these properties in GCE simulations. 
Specifically, as we have shown, one can apply the trans- 
formation of Eq. ([7]) to near-critical GCE data and then 
proceed as if the density distribution had been obtained 
within the RGE. The subsequent analysis in terms of the 
scaling of iso-q* maxima is arguably simpler to imple- 
ment than a widely used approach for GCE data [B] in 
which one determines a symmetrized ordering operator 
distribution via a linear field mixing approximation, and 
matches this distribution to a universal form to obtain 
effective critical parameters. 

Finally, with regard to the wider context of this work, 
it should be stressed that RGE simulations cannot be re- 
garded as competitive with the GCE in situations where 
the latter operates efficiently, e.g. for single component 
fluids or mixtures of similarly sized particles. Specifi- 
cally, the lack of a chemical potential field implies that 
histogram reweighting cannot be implemented to scan a 
range of overall densities or concentrations on the basis 
of a single simulation run; instead a separate RGE sim- 
ulation is required for each density po of interest. Where 
we expect the RGE to come into its own, is for highly 
size-asymmetric mixtures, for which GCE sampling be- 
comes ineffectual. As shown in ref. [T] the Generalized 
Cluster algorithm, when combined with the RGE, does 
permit efficient sampling of configuration space in such 
systems, and thus the methods we have described should 
facilitate the accurate determination of phase coexistence 
properties. We have already begun applying them in this 
context and hope to report on our progress in due course. 
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